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Abstract 

We present a new method for the numerical calculation of canonical reaction 
rate constants in complex molecular systems, which is based on a path integral 
formulation of the flux-flux correlation function. Central is the partitioning of 
the total system into a relevant part coupled to a dual bath. The latter consists 
of two parts: First, a set of strongly coupled harmonic modes, describing, for 
example, intramolecular degrees of freedom. They are treated on the basis of a 
reaction surface Hamiltonian approach. Second, a set of bath modes mimicking 
an unspecific environment modeled by means of a continuous spectral density. 
After deriving a set of general equations expressing the canonical rate constant 
in terms of appropriate influence functionals, several approximations are intro- 
duced to provide an efficient numerical implementation. Results for an initial 
application to the H-transfer in 6-Aminofulvene-l-aldimine are discussed. 
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INTRODUCTION 



Sophisticated experimental methods nowadays provide a rather detailed insight into 
molecular dynamics, unraveling the importance of quantum effects even in rather 
complex systems at room temperature^. This provides a challenge to theory since a 
fully quantum mechanical description of condensed phase dynamics remains out of 
reach and therefore approximate methods have to be developed. Among the oldest 
problem is that of chemical reaction rates which in fact gives a straightforward means 
for identifying quantum tunneling in terms of the non-Arrhenius behavior. Different 
methods have been developed to account for quantum effects in rate calculations (see, 
e.g., reviews in Refs. |2|31 which set the focus on enzyme reactions or the recent devel- 
opments in Ref. H]). A rigorous formulation of quantum mechanical rate constants can 
be given on the basis of the path integral approach'^', see also the related instanton 
type approaches, e.g. in Refs. I8HTU1 

Canonical rate constant are commonly calculated using the flux-flux correlation 
approach 5 , which requires a path integral propagation in complex time. Here, a break- 
through in numerical efficiency has been the quasi-adiabatic propagator (QUAPI) 
approach developed by Makri and coworkers^^. For the case of a generic system- 
bath model the QUAPI approach is based on a propagator splitting where the quasi- 
adiabatic path along which the bath oscillators are at their minimum position along 
the reaction path serves as the reference. In the context of rate calculations it has been 
applied to the situation of a double well, bilinearly coupled to a harmonic bath^^^l, 
and to electronically nonadiabatic reactions in Ref. [HI In another application Makri 
and Forsythe^ used the all Cartesian reaction surface Hamiltonian approach^ to 
determine a system-bath Hamiltonian for H diffusion in a silicon lattice. Employing 
a flexible bath reference for the Si environment led directly to the form of the Hamil- 
tonian used in the QUAPI method. However, to account for the two-dimensional 
motion of H an effective one-dimensional Hamiltonian had been used which was sup- 
plemented by an orthogonal harmonic mode with position-dependent frequency. The 
Si lattice bath modes were treated at the transition state geometry, i.e. mode-mode 
coupling and coordinate-dependence of the Hessian were neglected. In a subsequent 
publication, the issue of coupled bath modes has been addressed for a generic sys- 
tem's 1 . 

In the present work we consider the more general situation, where a large am- 
plitude reaction coordinate belongs to some polyatomic molecule, which is further 
embedded in some environment such as a solvent or a solid state matrix. The term 
polyatomic molecule is assumed to include situations with strongly coupled solvation 
shells. This setup will be termed system coupled to a dual bath. For the case men- 
tioned the distinction between intra- and intermolecular baths is motivated by the 
following observation: Quite often one faces a situation where the (intramolecular) 
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reaction coordinate is strongly coupled to specific intramolecular modes with an in- 
teraction potential that is not of the standard bilinear form. This coupling can be 
well-described by a reaction surface Hamiltonian, which in principle is amenable to an 
ab initio treatment. For the surrounding solvent this level of sophistication is often 
not necessary as the spectral densities associated with the coupling are broad and 
featureless. This suggests a treatment in terms of empirical models or classical calcu- 
lations of respective correlation functions, e.g., on the basis of molecular mechanics 
force fields^. 

The goal of the present paper is to develop a path integral expression for the 
calculation of canonical rate constants for a reaction coordinate coupled to a dual 
bath. This approach is applied to the case of H-transfer in 6-Aminofulvene-l-aldimine 
embedded in some model environment. Although our results for this case are of 
preliminary character, this reaction in principle shows some interesting effects. It 
was investigated in detail using NMR spectroscopy by Limbach and coworkers^. 
The temperature-dependent rate was found to be sensitive to the phase of the sur- 
rounding medium, which was either amorphous or crystalline. In general the reac- 
tion in the amorphous phase proceeds faster and the observed kinetic isotope effect 
(KIE) becomes temperature independent for low temperatures; at T =298 K it was 
k n /k B = 4. In contrast for the crystalline environment the KIE was temperature 
dependent throughout the measured range, which did not include tunneling regime; 
at T =298 K it was k R /k D = 9. The analysis of the experimental data was performed 
using the Bell-Limbach modeP^. This model introduces the reorganization energy 
for H-bond compression which is necessary for tunneling to occur from the intrinsic 
barrier for the transfer in the compressed state assuming a two step process. Further, 
a heavy atom mass effect is assumed for the transferred particle. Based on this model 
the effective barrier was estimated to be 3 kcal/mol and 1.9 kcal/mol for the crys- 
talline and amorphous phase, respectively. In both cases the reorganization energy 
amounted to 0.5 kcal/mol and the mass effect was found to be 1 a.m.u. Accounting 
for zero-point energy effects yielded an effective barrier for D transfer of 1.2 kcal/mol 
and 0.7 kcal/mol for the crystalline and amorphous phase, respectively. Although this 
model is of empirical character it shows the importance of the specific coupling to 
bond-compressing intramolecular modes as well as the influence of the environment 
on the reaction rates, thus illustrating the essence of the present dual bath approach. 

In the following we will start by introducing the system-bath Hamiltonian; a brief 
summary of the derivation of the intramolecular reaction surface Hamilton is given in 
the Appendix. Afterwards the path-integral expression of the canonical rate constant 
will be derived and some approximations simplifying the numerical treatment will be 
introduced. Subsequently, the application to the H/D-transfer in 6-Aminofulvene-l- 
aldimine is discussed and we conclude with a summary. 
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THEORY 



System-Dual Bath Hamiltonian 

Large amplitude motions of certain coordinates of a polyatomic molecule embedded 
in some environment will be described as a relevant low- dimensional system coordi- 
nate, s, coupled to a dual bath, where the intramolecular and environmental degrees 
of freedom (DOFs) are denoted Q and q, respectively. In the Appendix we give a 
brief account on the derivation of a reaction surface Hamiltonian for the intramolec- 
ular problem,^ specified to the case of a linear reaction path^. The resulting in- 
tramolecular Hamiltonian can be written as the sum of the (one-dimensional) reaction 
coordinate part (we use mass-weighted coordinates and atomic units throughout) 

H o--l^ + Vo(s), (1) 

an intramolecular bath part 

*=iE[-|s+<*a]. p) 

and a coupling part 

Vi( 8 , Q) = -j2ft (s) Qk+lYl Kkki ^ - \ E ^ ( 3 ) 

k k,k' k 



Here, f(s) is the vector of forces exerted on the oscillators (Eq. (A. 10)), K(s) is 
the reaction coordinate dependent force constant matrix, and io\ = K kk (s Te {) is the 
frequency of the kth bath mode at some reference value of the reaction coordinate. 

The coupling of the reaction coordinate and the intramolecular DOFs to the har- 
monic bath of the environment, 



"» = sE[-|> +<■*«:]■ (4) 



a ' a 



will be assumed to include the lowest-order terms of a Taylor expansion with respect 
to q, i.e. 

{s)q a + S ^C^ k (s)Q k q a . (5) 



a,k 



Here, where d a (s) and c a< k(s) are some coupling functions to be specified for the 
system at hand. Thus the total Hamiltonian is given as 

H = H (s) + H^Q) + H 2 (q) + V x (s, Q) + V 2 (s, Q, q) . (6) 
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Canonical Quantum Reaction Rate 

We will use the flux-flux correlation function expression of the reaction rate between 
reactant and product, k^p, due to Miller and coworkers^ 

kRP = ^J™C { (t)dt, (7) 

where 

C t (t) = Tr {Fe iHt 'Fe~ lHtc } (8) 

is the autocorrelation function of the symmetrized flux operator specified here to the 
case of a one- dimensional reaction coordinate s with the dividing surface at s = 0, 
F — | (p s 8(s) +S(s)p s ), and Z = Tr(— {3Hn) is the canonical partition function of 
suitably defined reactant Hamiltonian Hr. The complex time, t c = t — i(3/2, is due 
to the combination of the time evolution operator and the Boltzmann operator and 
P = 1/kaT. 

The flux autocorrelation function can be calculated by approximating the momen- 
tum operator in the vicinity of the dividing surface by a finite difference expression 
with increment As^ 

°t (*) = l K ( As > As > °> °> *c) - K{0, As, 0, As, t c )\ , (9) 

where 

/oo 
dQdq(Q\(q\(s' ,/ \e iHt ^s / ')(s / \e-^\s)\q)\Q). (10) 
■oo 

The elementary propagators in this expression can be evaluated using the path inte- 
gral technique, i.e. dividing the complex time t c into N slices. This yields^ 

K(si, SjV+l, SjV+2, S2N+2, t c ) 

dQdqds2 ■ ■ ■ dsNdsN+3 • ■ • ds2N+i 



oo 

N+2 



x(Q\(q\ JJ (s n+1 \e- tH6 "\s n }]J(s n+1 \e- tHS -\s n }\q}\Q} 

n=27V+l n=N 

/oo /*oo 
■ ■ ■ I ds2 ■ ■ ■ dSNdSN+3 ■ ■ ■ d,S2N+lFmfl(si, «2, ' " " , S2N+2, tc) 
-00 J — 00 

N+2 1 

J] (s n+1 |e"^^|s n ) n(s n+1 |e-^^|s n ), (ll) 



x 

n=2iV+l n=iV 
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where the time steps S n are defined as follows: 

-t: 



&2N+2 

5„ 



JN+2 



2N 



—t* 
N 



61 

to 

TV' 



, n = jV + 3,- 

tc_ 

2N 

, n = 2, • • • ,N. 



2N + 1 



(12) 



In Eq. (11) the influence functional is defined as 

/oo [ 
dqdQ(q\(Q\ j ] e - i [^( < 3)+^( s ».Q)+^(9)+^(^ > Q,gFn|Q)|q) ) 
"°° n=2N+2 

(13) 

where {s n } denotes a specified path realization. With the help of the exact prop- 
agator for harmonic oscillators^ one can obtain the following result, e.g., for the 
intramolecular bath part 



[Q n+1 \e-< H ^ +v ^ Q ^\Q n ) = expH^QJUl! 



2iri sm(u> k 5 n ) 



x exp 



E 



cot(u: k S n )Q nk 



> Qn+l,kQnk 



sm(u k S r , 



■ (14) 



Note that Vi(s n ,Q n ) not only contains the force on the oscillator coordinates but 
also the mode-mode coupling and the change of the diagonal elements of the force 
constant matrix with respect to the chosen reference value of the reaction coordinate. 
Actually the choice of the latter does not play an important role, if one assumes that 
5 n is chosen to be sufficiently small. 

Using a similar expression for the environmental part of the Hamiltonian, we arrive 
at the following influence functional 

/oo 
dQ x . . . dQ 2N+2 dq 1 . . . dq 2N+2 exp{g({s n }, Q, q)} (15) 
■00 



g({s n },Q,q) 



E 

nk 



lUJ k 



cot(uj k 5 n ) + ^ fc r f'" y ) Q': nk 



a ) \ / v2 Qn+l,kQnk 



sm(u k S n ) 



+i ^ 5nfk{Sn)Qnk ~ ~ ^ QnkK kk '(s n )Q nk > 



nk 



ikk' 



' sin w„o, ■ L 



nka 



(16) 



where 

Fq = 11 J 2msm(uJ n ) (17) 

and 



nfc ' 

are path independent prefactors. 

The partition function can be calculated following the same lines 

/oo /»oo 
• • • / ds 1 ds 2 ■ ■ ■ ds Np F )3 (s 1 , s 2 , ■ ■ ■ , s N/3 ) 
-oo J — oo 



-oo J — oo 

1 



x ( Sl \e- iH °^\s N0 ) J] (s n+1 \e- iH °^\s n ) (19) 

n=N p -l 



with the influence functional 

i 



/OO -i 
dqdQ(q\(Q\ ]^ e - i[Hl(g)+yi(s "' Q)+H2(g)+y2(s "' Q ' q)] ^|Q)|g). (20) 

Here 5p = —ift/Np and Np = 2N is the number of time slices for the imaginary time 
-ifi. 

In a next step we need to evaluate the integrals in Eq. (15) which are of the 
following complex-coefficient Gaussian type 

/oo /»oo 
••• / dx 1 dx 2 ---dx N eyv{-^A mn x m x n + i^B mn x m x n + ^W n x n }, (21) 

where both A and B are real symmetric matrices. For any physically meaningful 
case the matrix A is positive-definite and the integration will converge. Moreover, it 
is possible to find one invertible real matrix, U c , to congruently diagonalize A and B 
simultaneously. If Ui and U2 are orthogonal matrices such that 

Ui T AUi = a = diagjai, a 2 , ■ ■ ■ , a^} 
U 2 T a^Ui T BUia^U 2 = b = diag{&i, b 2 , ■ ■ ■ , b N } (22) 

and the matrix A is positive-definite such that all eigenvalues {a n } are positive, one 
can define the transformation matrix 

U c = Uia^u 2 (23) 

which transforms UjAU c = 1 and UjBU c = b. Using the new variables {y n } 
defined by y n = ^^(U" 1 )^^ the integration in Eq. (21) can be performed analyt- 
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ically to give 

/OO POO 
■■■ dx 1 dx 2 ■ ■ ■ dx N exp{- ^ A mn x m x n + i ^ B mn x m x n + ^ W n x n } 
00 J -00 mn mn n 

/oo /»oo 
• • • / c^cft/a -~dy N exp{- - i& n )0* + ^ w nVn} 

(24) 



n 




1 



1 — ib, 



cxp 



4(1 - ib n ) 



where w n = ^ m (Uc)nmWm. Here and in the following the square root of a complex 
number means its principal value, i.e., the non-negative real part. 

Using this method it is at least in principle possible to solve Eq. (15). However, in 
practice this would imply to numerically diagonalize a large matrix for each specified 
path. In order to simplify matters we reconsider the environmental bath part. Here, 
the quadratic coefficients of the bath oscillators are path independent and assumed 
to be uncorrelated between each other. Based on above mentioned procedure we can 
find a frequency-dependent real invertible matrix U g (u) to congruently diagonalize 
each bath mode 

Qn = y^J^g 1 ((*i)]nn'Qn' (25) 



such that 



V] ■ - f * x [cos(u5 n )ql - q n+1 q n ] = - V(l - ib q n {uj))q 



n ! 



(26) 



where the {b^iuj)} (and (a^(o;)} which will appear below) are the eigenvalues from di- 
agonalizing the corresponding coefficients matrix according the procedure introduced 
in Eq. (22). Using the new variables q na = Y^n'^q 1 ( u a))nn'Qn'a the integration over 
{Qna} can be performed analytically. The final result for influence functional is given 
by 



Fm ({>«}) 

g({s n },Q) 



F Fr,F 



dQ 1 dQ 2 ■ ■ ■ dQ N exp{g({s n }, Q)} 



nk 



lbJ k 



cot(a; fc 5 n ) 



Qnk 



Qn+l,kQ 



nk 



sin(o;fc5 n ) _ 



+i^2<> n f k (s n )Q nk - ^ ^ QnkK k k'(Sn)Qnk' + A({s n }) 
nk nkk' 

Qnk,n' k 1 [Sni Sn'^QnkQn'k 1 •> 

nk nkn'k' 



where 



n 



7T 



a n (u} a ) V (1 - ib n (uJa)) 



(27) 



(2f 
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A(K}) 



2-^ A\i _ ,7 



^n'a,nk \^n) 
9nk,n' k 1 i^Sni ^n') 



4[1 - z&£(w a )] 
/J[U g (a; Q ,)] n / n 5 n 'rf a (s n /) 

n' 

W n ' a U n ' a n k 

2[1 -ib%(u a )] 



-i[U q(u Q )] nn >5 n C ak (s n ) 



E 



E 



U>n" a,nkV>n" a,n' k' 



4[1 



(29) 



So far the g-integrations have been performed following the idea from Eq. (21) 
to Eq. (24). The quantities entering Eq. (27) can be calculated readily before the 
Q-integrations. In a final step the procedure of Eq. ( 22 ) can be applied to the Q- 
integrations to numerically diagonalize the complex coefficient matrix in Eq. (27) for 
each path of the reaction coordinate. The final result for the influence functional of 
the reaction coordinate plus dual bath system can then be formally written as: 



Finfl (K» = F q F Q F q e A ^ H 



nk 




cxp 




(30) 



with the different functions defined in Eqs.(17), (18) and (28). The quantities a® k , 
b® k , and w nk in above formal expression can be obtained from the numerical diago- 
nalization of the respective complex coefficient matrix. 



Approximations 



Depending on the system size obtaining the quantities in Eq. (30) by direct diagonal- 
ization might become rather time consuming, due to those terms which depend on 
the system's coordinate and, therefore, have to be evaluated for each specific path. 
Therefore, we will introduce certain approximations to make the approach numerical 
efficient for such cases. 

First, we will assume that the coupling strength between the Qk and q a modes does 
not strongly depend on s. Thus we ignore the s-dependence of the coupling strength 
between Q k and q a , i.e, {C Q k} are simply constants and hence {g n k,n'k'( s n, s n') — 
g n k,n'k'} are also constants. Next, we assume that not for all modes, {Qk}, the mode 
mixing due to K(s) shows a strong coordinate dependence. In the following we 
use {Q k } to denote those intramolecular DOFs, which are most strongly affected by 
the motion of the reaction coordinate s. The remaining intramolecular modes are 
comprised in {Q u }- For the latter modes, the quadratic coefficients will be replaced 



by their s-independent mean values along the reaction path in Eq. (27), i.e., 



1 



L 



K vv ,{s n ) {K vv ,) = — / K w ,(s)ds, (31) 
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where 2L is the length of the reaction path. Under this approximation, we need to 
diagonalize a large matrix just once, while for each specified path we only need to 
diagonalize a much smaller matrix since only a few DOFs, {Qk}, are significantly 
coupled via s. 

Following the idea of Eq. (22) we can find a real invertible matrix which 
congruently diagonalizes the quadratic coefficient matrix related only to {Q v } 



Ldyb~n) \ ^2 Qn+l,uQr 



sm(co u S r , 



'^iujy \cot(u v 5 n ) + v Y^ )Q 

nu L ^ 

— n 5^ ^ n i-K-vv' ) QnuQnu' + 9nu,n'u'QnuQn' , 



2 

nui 



(32) 



where Q nu = Yln'v' nun'u' Q nv> ' With the help of this transformation we can 
analytically integrate over the {Q n u} part. This will further contribute a pre-factor 
Fq and some modifications to the exponential factor compared with Eq. (27) 



-Finfl {{Sn}) 
g({Sn,Qnk}) 



F q F Q F q F Q 



dQ 1 dQ 2 ■ ■ ■ dQ N exp{g({s n , Q nk })} 



nk 



%UJ k 



<zoX>(oj k b n ) + 



(u k 5 r , 



Qnk 



Qn+l,kQ 



nk 



sm(u k 5 n ) 



+i ^ $nfk{s n )Qnk ~ ^ ^ QnkK k k> {s n )Qnk> + A({s n }) 
nk nkk' 

+i 2_j $n^fk{Sn)Qnk + 9nk,n' k'QnkQn'k' + A({s n }) 

nfc nkn'k' 

+i ^ Sn^fk{Sn)Qnk + 9 nk,n' k'QnkQn'k' j 

rafc nkn'k' 



(33) 



where 



n 



7T 



',,,/ y 1 ib nv 



(34) 
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and the additional terms caused by the reduction of DOFs are defined as follows: 



A(K» 
w nu 
A/ fc (s n ) 

Unk,n'u' 
Qnk,n'k'\^m ^n'J 



^ 4(1 - ife^) 

nV ,nu 



E UJn> 
, , 2(1 - i6 nV ) 

n"u,n'u' Qnk,n"u 1> 

5 n K ku (s n ) 



fJ"nk,n"u^"n'k',n"u 



t, 4(1 -zVi/) 



Similar to Eq. (30) the final result can be written formally as 



Ftea ({*„}) = F g F Q F g F Q e A «^» +A « s "» J] 




«nfc V 1 ~ i6 nfc 



exp 



(35) 




with the different functions defined in Eqs.(17), (18), (28), and (34). The quantities 



Vfc> nki 



and w n k in above formal expression can be obtained from the numerical 
diagonalization of the complex coefficient matrix and the final sum is only for modes 



which strongly couple to s. The final numerical calculations may start from Eq. (33) 



which is feasible since only a very low-dimensional matrix (according to the coordi- 
nates {Qk}) needs to be diagonalized for each specified path. 



APPLICATION TO THE H/D-TRANSFER IN 6-AMI- 
NOFULVENE-1-ALDIMINE 

In this section we present results of a preliminary simulation based on a reaction sur- 
face model Hamiltonian describing the intramolecular H atom transfer in 6-Aminofulvene- 
1-aldimine. Here, our aim is not to provide a quantitative assessment of this reaction, 
but to illustrate the theoretical formalism presented in the previous section. The 
configuration of two stationary points, which have been obtained at the B3LYP/6- 
31+G(d,p) level of theory^ 1 are shown in Fig.[l]. The minimum configuration in panel 
(a) corresponds to the reactant or equivalent product. The H atom transfer process 
can take place from the reactant via the transition state (panel (b)) to the product 
or inversely. The reaction barrier height, as calculated by the energy difference of the 
minimum and the transition state, is 3.8 kcal/mol (fully relaxed gas phase barrier). 

The unit vector which defines the linear reaction path is given by the direction 
pointing from the equivalent reactant to the product, i.e., e s = (i? pro d— -R rC ac) / (l-Rprod— 
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Figure 1: 6-Aminofulvene-l-aldimine: (a) minimum configuration, (b) transition state 
for the hydrogen atom transfer, (c) Strongly coupled normal mode at the transition 
state having a frequency of 641 cm -1 . All results have been obtained at the B3LYP/6- 
31+G(d,p) level of theory. 

-Rread) (cf. Fig. [I]). The potential along the one-dimensional linear reaction path co- 
ordinate, V(Rq), is shown in Fig. [2] According to the present linear reaction path, 
the barrier is as high as 14.85 kcal/mol. 

Within the harmonic approximation for the intramolecular bath modes, the energy 
difference with respect to the fully relaxed reaction path will be recovered by the so- 
called bath reorganization energy^ (see also discussion in Refs. 123 and [23 where 
two and three reaction coordinates, respectively, have been used to obtain a better 
agreement with the fully relaxed barrier even without taking into account coupled 
harmonic vibrations). Furthermore, in the real system, there will be a contribution 
to the reorganization energy due to the interaction with the environment. In the 
following we do not attempt to fit the environmental contribution such as to obtain 
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-9 -6 -3 3 6 9 
Reaction coordinate s (As) 

Figure 2: The zeroth-order potential energy curve, V(Rq), obtained at the B3LYP/6- 
31+G(d,p) level of theory for the hydrogen/deuterium atom transfer reaction in 6- 
Aminofulvene-l-aldimine as shown in Fig.[T] The curve has been generated from 19 
points, symmetrically distributed with respect to the reference geometry s re f = 0, 
which corresponds to i? re f = (-R pro d _ -Rrcac)/2. The react ant /pro duct configuration 
is at s = ±6As, i.e. As = |-R pro d — -Rread/12. 



agreement with the experimental estimate by Limbach and coworkers.^ 

The considered molecule has 105 intramolecular vibrational degrees of freedom, 
whose couplings to the reaction coordinate, i.e. fk(s), and Kk,k'{s), can be obtained 



as described in the Appendix, Eq. (A. 18). For the reference i£ re f we have chosen the 
point midway between reactant and product along the one-dimensional reaction path. 
For the present illustration we have selected only one strongly coupled mode Qk for 
explicit consideration. The displacement vectors are shown in Fig. ^p. Apparently 
this mode symmetrically modifies the H-bond length and therefore modulates the re- 
action barrier. The remaining intramolecular modes as well as possible environmental 
modes are comprised into the bath q. In other words, we have simplified matters and 



started directly from Eq. (27) 



The coupling between the environment q and the intra-molecular DOFs, s and 
Qk, are defined as 

d a (s) = d 1 e~^ ( %(s + 7] s 2 ) 
C a , k (s) = C a , k = c ie -^-^ 2 ^, (36) 

where d%, d 2 , co , T], ci, and c 2 are parameters. The s-dependence of d a (s) has been 
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4 l i i i i i i 

2.5 3 3.5 4 4.5 5 5.5 6 

1000/T (K _1 ) 

Figure 3: The calculated temperature dependence of H/D transfer rate constants, 
Eq. ([7]), in the thermal activation region based on a one- dimensional linear reaction 
path, s, coupled to one intra-molecular mode, Qk, and 50 bath modes, q. In the 
simulation the reaction coordinate integration has been replaced by a sum over the 
three points, s = —1, 0, 1. For the partition function of the reactant these points were 
chosen as s = —7, —6, —5. The number of time slices has been N — 4. 

expanded to second-order and the bath frequency dependence has been simply chosen 
to be of Gaussian form. The environmental modes are assumed to have uniform 
density of states in the region, where we take into account the coupling with the 
molecular DOFs. 

The calculated canonical rates, Eq. ([7]), obtained from this preliminary model 
Hamiltonian are shown in Fig.[3] At 298 K the KIE is kp P /kp P = 10, when the 
following coupling parameters are used: c\ = di =(0.628 kcal/mol) 2 , C2 = c?2 = 6.28 
kcal/mol, and rj = 0.2As _1 . The involved bath frequency region covers the range 
from 3 to 30 kcal/mol with 50 harmonic oscillators equally distributed. Given the fact 
that the experimental KIE ranges between 4 and 9 and strongly depends on the phase 
of the environment, the present order-of-magnitude agreement is rather reasonable 
given the simple model for the system-bath coupling. Further, we note that the same 
holds true for the absolute values of the rates. However, the obtained values for the 
activation energies (slopes of curves in Fig. [3]) deviate from the experimental ones. In 
Ref. [19] it was found that the rate between the thermal activation energies between 
the H and the D case is about 2/3. In the present simulation the activation energies 
are about 3-4 times the experimental ones and the difference between the activation 
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energies of different isotopomers are too small. The latter fact is not surprising since 
the shapes of potential curves for hydrogen and deuterium transfers are the same 
and the only difference lies in the length of the step As which appears in Eq. Q. 
The ratio for the steps is only slightly different from one, As(H)/As(D) = 0.9978. In 
order to improve the description at this point, more intramolecular vibrational modes 
need to be taken into account. This would lower the effective reaction barrier due 
to reorganization energy contributions and therefore the activation energy. On the 
other hand, due to the isotope dependent effective coupling, the difference between 
H and D activation energies would become more pronounced. 




2.5 3 3.5 4 4.5 5 5.5 6 2.5 3 3.5 4 4.5 5 5.5 6 2.5 3 3.5 4 4.5 5 5.5 6 
1000/T (K" 1 ) 1000/T (K" 1 ) 1000/T (K _1 ) 



Figure 4: The left and middle panels show the convergence of the thermal activation 
energy (the slope) for H-transfer in the high temperature region by only considering 
configurations which are important for thermal activation, i.e., around s = 0. In the 
left panel the dependence on the number of time slices N is shown for s = —1, 0, 1 
and in the middle panel the case s = —2, —1, 0, 1, 2 is given for N = 4. The right panel 
shows the quantum tunneling effects in the low-temperature region by covering some 
configurations which are important for tunneling (s = —2, —1,0, 1,2 (dashed curve) 
and s = —5, —1, 0, 1, 5 (solid curve)). 



The numerical effort in the calculation of the propagator in Eq. (11) by path in- 
tegration depends on the number of time slices, N, as well as on the method for 
evaluating the multi-dimensional integrals of the reaction path coordinates, s n . For 
the present application in Fig.[3]the focus has been on the thermal activation range, 
i.e. the high-temperature regime. This allowed us to simplify the rate calculation 
by performing the integration as a sum over three points in the vicinity of the reac- 
tion barrier, s = —1,0, 1. In Fig. [4] (middle panel) we show the dependence on the 
number of discretization points for N = 4 time slices. Specifically, we have chosen 
s = —2, —1, 0, 1, 2. The ignorable difference shows the applicability of the simplifica- 
tion technique, which we have adopted for the high temperature calculations. The 
dependence on the number of time-slices is shown in the left panel of Fig. |4| Clearly, 
the variation for the covered range, N = 2,3,4, is rather small, justifying our choice 
of iV = 4 in Fig. [3] Finally, we address the issue of tunneling in the right panel of Fig. 
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|4| In principle, accounting for quantum tunneling at low temperatures requires to 
include configurations, which are located near the turning point corresponding to the 
energy of the tunneling particle. In order to illustrate this point, we present results for 
five discretization points, i.e. s = —2, —1, 0, 1, 2 (dashed curve) and s = —5, — 1, 0, 1, 5 
(solid curve). The change of the mechanism from thermal activation to tunneling is 
apparent from the quite different temperature dependence of the rates. In passing we 
note that a systematic study of different discretizations might give an indication for 
those configurations that contribute to the tunneling process. 

SUMMARY 

We have developed a path integral method for the determination of canonical reac- 
tion rates for the case of a reaction coordinate coupled to a dual bath. The latter 
is comprised of an intramolecular part, which is modeled using the reaction surface 
Hamiltonian approach, and an intermolecular (solvent) part. Such a partition of the 
interactions into different structurally motivated levels appears to be most suitable 
for the description of intramolecular proton or H-atom transfer reactions. The formu- 
lation benefits from the harmonic oscillator nature of the intra- and intermolecular 
baths in two respects: First, it enables us to perform the integration over the bath 
variables and second the reaction surface Hamiltonian method provides a means to 
determine intramolecular Hamiltonian parameters from first principles. This involves 
couplings between normal modes along the reaction path due to the non-diagonal 
Hessian matrix, which require a diagonalization for each specified path and thus sub- 
stantial numerical effort. We have suggested an approximation which amounts to the 
replacement of the reaction coordinate dependent Hessian by its averaged value for 
less strongly coupled intramolecular modes. 

The initial application has been to the H/D transfer in 6-Aminofulvene-l-aldimine. 
Despite the various additional approximations in this application it could be shown, 
that our approach can give reasonable (i.e. order-of-magnitude) estimates for the 
reaction rates. Further work on this particular reaction shall be directed to obtain 
an improved description of the reaction surface as well as a realistic model for the 
interaction with the environment. 
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APPENDIX 



Reaction Surface Hamiltonian 

The reaction surface Hamiltonian combines the description of several large amplitude 
coordinates {s a } coupled to many small amplitude displacements j ESEH28I rp Q 
generate this Hamiltonian from the exact Cartesian coordinate Hamiltonian we can 
directly exploit our recently developed kinetic energy quantization method.^ 
Suppose that we have the Cartesian Hamiltonian 

H(R) = T(R) + V(R) 



1 1 d 2 

T(R) = -P 2 = %, (A.l) 

y ! 2 2 OR 2 K J 

where R is the 3iV-dimensional vector of mass-weighted Cartesian coordinates for 
system with N atoms and P = —id/dR is the corresponding linear momentum 
operator. Suppose that there is a reaction surface defined by a function along the 
reaction coordinates s, i.e., 

R = R (s) . (A.2) 

The potential energy function, V (R), is expanded around the reaction surface as 
follows 



dV 



1 . ^d 2 V 



V(R)^V(R (s)) + AR(s) L — + -AR(s) L ^ AR(s) , (A.3) 



n ... 2 OR 1 



R Q (S) 



where AR (s) = R — R (s). The reaction surface is defined in such a way that the 
potential energy V (R) can be approximated by low-order orthogonal displacements, 



i.e., Eq. (A.3) can be truncated in the given form. 

To obtain the reaction surface Hamiltonian we first need to define the new coor- 
dinates, i.e., the reaction coordinates {s a } and the orthogonal displacements {Qk}- 
The former are already defined by the reaction surface as well as the unit vectors 
{e a (s)} according to which we have the reaction coordinate vector 

D 

s = J2 s * e * ( s ) • ( A - 4 ) 

a=l 

To get the latter we need a projection operator to project out the reaction coordinate 
s 

V(s)^l-J2^e T a . (A.5) 

a 

Then we can diagonalize the projected Hessian matrix K (s) for each point of the 
reaction surface by an orthogonal transformation Urs (s) 

XJ RS (s) f K (s) U RS (s) = diag{- • • a£ (a) ■ ■ ■ u 2 (s) ■ ■ ■ u 2 (s) ■ ■ ■ }, (A.6) 
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where K (s) = V {s)d 2 V/dR 2 \ R V (s) is a real symmetric matrix. 

In total there are D + 6 zero eigenvalues {u^} and {oj 2 } corresponding to the 
reaction coordinates and six-dimensional global translation and rotation, respectively. 
The orthogonal transformation matrix contains the corresponding eigenvectors of 
K(s) 

U RS (s) = (• • • e a (s) • • • e g (s) • • • e k (s) • • - ) . (A.7) 

The six-dimensional global translation and rotation as well as the 3N — 6 — D dis- 
placements orthogonal to the reaction surface are defined by 

R g = e]AR 

Q k = ejAR. (A.8) 
The original 3Af-dimensional vector is now expressed with the new unit vectors 

R = R rel + ^2s a e a + ^2R g e g + ^2Q k e k , (A.9) 

a g k 

where the reference geometry R Te { = R (s = 0) is the origin of the new coordinates 
system. 

Based on the knowledge of the new coordinates it is not difficult to find the 
potential energy 

V (s, Q) = V(R {s))-J2 fk (a) Q k + lj2 Uk ^ (A.10) 

k k 

where f k (s) = —eJdV/dR\j^. It is obvious that the potential energy does not 
depend on {R g }, however, the kinetic energy operator (KEO) does depend on {R g } 
and normally it is not possible to separate them exactly. According to Ref. [29] the 
following formal KEO can be obtained 

m l=jdR f dR\ T ~ , A . 

T= 2 P m[M) p ' (A ' U) 

where R = ( s T R g T Q T ) is the full set of the new coordinates and P = 
—id/dR. According to Ref. [29] all components of P are Hermitian due to the or- 



thogonality of transformation except P§. Eq. ( A.ll ) has a fully coupled form in case 
of a general reaction surface. The factor which is responsible for complexity when it 
comes to a numerical implementation is that all unit vectors depend on s, i.e., the 
orthogonal transformation matrix Urs (s) depends on s thus we have to calculate the 
derivatives with respect to s. 
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Linear Reaction Surface Hamiltonian 



In the following we will simplify the KEO, Eq. (A. 11), by choosing a different repre- 
sentation in terms of constant unit vectors that describe the reaction coordinates s 
for the special case of a linear reaction surface.- 22 ' With the help of certain predefined 
constant unit vectors {e a } we can obtain the following equation for the linear reaction 
surface 

R (s) = R Tei + s * e a . (A. 12) 

a 

The coordinate transformations are modified as follows: 

R = R Q (s) + Q k e k = R rei + s a e a + Q k e k 



elAR, Q k = elAR, 



(A.13) 



where AR = R — R re f is different from AR (s) in Eq. (A. 3 ) while {Q k } and {e k (s)} 
have the same definition as in the previous section. Note that we have combined the 
{R g } and {Q k } into the same set of indexes {Q k } to simplify the notation. With 
the help of Eq. (A. 11) and Eq. (A.13) we can derive a simplified KEO for a linear 
reaction surface. First, we calculate the elements of the Jacobi matrices starting from 
Eq. (A.13). Using the chain rule to calculate the derivatives from Eq. (A.13) leads to 

I = e a 

>de k s 



the following results 



dR 

dQ k 



dR 



AR 1 



dSr 



(A. 14) 



Thus the elements for the matrix products in Eq. (|A.11|) can be obtained as follows 
ds ( <9s N ' 



- 



dR \dR 
ds fdQ 



Sat 



dR \dR 

dQ (dQ_ 
dR [dR 



a/3 



ak 



kk' 



'kk' 



£ (AR 




AR 



T de k 
ds n 



(A.15) 



Based on above equations we can simplify Eq. (A.ll) to yield (cf. Ref. |22 



T( 8 ,Q) 



2 y^ P a + 2 y^ P k ( Skk' + ^2 B <*k B ak' J Pk> 
a kk' \ a / 

+ ( \ Pa BakPk + h - c - ) ' 

V ak J 



(A.16) 
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where B ak = AR T de k /ds a . Note here all the components of momentum are Hermi- 
tian according to Ref. |29j The kinetic couplings are caused by the s-dependence of 
{e^} as can be seen from the expression for B ak . The potential energy is still given 



by Eq. flA.10[ ). 

The KEO can be further simplified by using more constant unit vectors for the 
expansion of the coordinate space, i.e., we get rid of the s-dependence of {efc} (for 
alternative approaches see also Refs. [22|29|. The most simple case, in which the 
kinetic energy has a quite trivial form while the potential energy is no longer diagonal, 
is the space whose unit vectors are all constants. This can be achieved by diagonalizing 
the projected Hessian matrix at only one point, -R re f, instead of each point on the 
reaction surface. The new representation is obtained by a pure s independent rotation 
and the new variables are defined by 

S a = €-a (R — -Rref) 

Q k = eJ(R-R rci ). (A. 17) 

Here {Qk} denote the remaining 3N — D variables which are the global translation, 
rotation and normal modes only at the reference point. Notice, that within this 
approximation overall rotations are not strictly projected out for a general point on 
the potential energy surface. The Hamiltonian in terms of the new coordinates reads 



2^ a 2 

a k 



I V — - - V — 



V(s,Q) = V{R Q )-Y J fk{s)Q k + \Y J K k k'{s)Q k Q kl , (A.18) 



2 

k k,k' 



where f k has the same definition as before and 

,d 2 V 



K k k' is) 



,T 1 



OR 2 



e w . (A.19) 

Jxo 



This form of the Hamiltonian has been used in the present paper to model the coupling 
between the reaction coordinate and the intramolecular vibrational modes in the 
application to the proton transfer in 6-Aminofulvene-l-aldimine. 
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